• We will use functions from packages base, utils, and stats (pre-installed and pre-loaded)
• We may also use the packages below (specifying package::function for clarity).

# Load pckgs for this R session
options(scipen = 999)
# --- General 
library(here)     # tools find your project's files, based on working directory
library(dplyr)    # A Grammar of Data Manipulation
library(skimr)    # Compact and Flexible Summaries of Data
library(magrittr) # A Forward-Pipe Operator for R 
library(readr)    # A Forward-Pipe Operator for R 
library(tidyr)    # Tidy Messy Data
library(kableExtra) # Construct Complex Table with 'kable' and Pipe Syntax

# ---Plotting & data visualization
library(ggplot2)      # Create Elegant Data Visualisations Using the Grammar of Graphics
library(ggfortify)     # Data Visualization Tools for Statistical Analysis Results
library(scatterplot3d) # 3D Scatter Plot

# --- Statistics
library(MASS)       # Support Functions and Datasets for Venables and Ripley's MASS
library(factoextra) # Extract and Visualize the Results of Multivariate Data Analyses
library(FactoMineR) # Multivariate Exploratory Data Analysis and Data Mining
library(rstatix)    # Pipe-Friendly Framework for Basic Statistical Tests
library(car)        # Companion to Applied Regression
library(ROCR)       # Visualizing the Performance of Scoring Classifiers

# --- Tidymodels (meta package)
library(rsample)    # General Resampling Infrastructure  
library(broom)      # Convert Statistical Objects into Tidy Tibbles

• Logistic regression is a classification model used with a binary response variable, e.g.:
  • yes|no, or 0|1, or True|False in a survey question;
  • success|failure in a clinical trial experiment;
  • benign|malignant in a biopsy experiment.
• Logistic regression is a type of Generalized Linear Model (GLM): a more flexible version of linear regression that can work also for categorical response variables or count data (e.g. poisson regression).
• When logistic regression fits the coefficients \(\beta_0\), \(\beta_1\), …, \(\beta_k\) to the data, it minimizes errors using the Maximum Likelihood Estimation method (as opposed to linear regression’s Least Squares Error Estimation method).

If we have predictor variables like \(x_{1,i}\), \(x_{2,i}\), …, \(x_{k,i}\) and a binary response variable \(y_i\) (where \(y_i = 0\) or \(y_i = 1\)), we need a “special” function to transform the expected value of the response variable into the \([0,1]\) outcome we’re trying to predict.

The logit function (of the GLM family) helps us determine the coefficients \(\beta_0\), \(\beta_1\), …, \(\beta_k\) that best fit this sort of data and it is defined as: \[ logit (p_i) = \ln\left( \frac{p_i}{1-p_i} \right)= \beta_0 + \beta_1 x_{1,i} + \cdots + \beta_k x_{k,i} \] where \(p_i\) is the probability that \(y_i = 1\)

Via the inverse-logit function (logistic), we solve for the probability \(p_i\), given values of the predictor variables, like so:

\[ P(y_i = 1 | x_{1,i}, \ldots, x_{k,i} ) = \frac{1}{1 + e^{-(\beta_0 + \beta_1 x_{1,i} + \cdots + \beta_k x_{k,i})}} \]

Then we examine the relationship between the binary outcome variable heart_disease and the binary predictor variable coffee_drinker.

• we use the handy count function from dplyr to count occurrences in categorical variable(s) combinations.

# Dataset manipulation
heart_train %>% 
  # count the unique values per each group from 2 categorical variables' combinations
  dplyr::count(heart_disease, coffee_drinker, name = "count_by_group") %>% 
  dplyr::group_by(coffee_drinker) %>% 
  dplyr::mutate(
    total_coffee_class = sum(count_by_group),
    proportion = count_by_group / total_coffee_class) %>% 
  dplyr::ungroup() %>% 
  # filter only those with heart disease
  dplyr::filter(heart_disease == "Yes_HD") %>% 
# Plot  
  ggplot(aes(x = coffee_drinker, y = proportion, fill = coffee_drinker)) + 
  geom_bar(stat = "identity") + 
  scale_fill_manual(values = c("#57b7ff", "#e07689")) + 
  theme_minimal() +
  ylab("Percent with Heart Disease") +
  xlab("Coffee Drinker (Y/N)") +
  ggtitle("Figure 3: Percent of Coffee Drinkers with Heart Disease") +
  labs(fill = "Coffee Drinker") + 
  scale_y_continuous(labels = scales::percent)

• In principle, we could use a linear regression model to study likelihood of having Heart Disease in relation to risk factors: \[ Y = \beta_0 + \beta_1 \cdot \text{(US\$ Spent Fast Food)} + \beta_2 \cdot \text{(Coffee Drinker = YES)} \]

• But we’ll see why the logistic regression model is a better option.

\[ \log \left( \frac{p}{1-p} \right) = \beta_0 + \beta_1 \cdot \text{(US\$ Spent Fast Food)} + \beta_2 \cdot \text{(Coffee Drinker = YES)} \]

• Let’s compare the 2 models (for simplicity, we ignore the coffee_drinker variable for now.):

# --- 1) Linear regression model
linear_mod <- lm(heart_disease ~ fast_food_spend# + coffee_drinker
                 , data = heart_data)
# --- 2) Logistic regression model
logit_mod <- glm(heart_disease ~ fast_food_spend# + coffee_drinker
                 , data = heart_data, family = "binomial")

We can now extract the coefficients from the linear regression and logistic models, then estimate the predicted outcomes from these models: lin_pred, logit_pred, and logistic_pred (i.e. the conversion form log(odds) to probability).

# --- 1) Extract coefficients from linear regression model
intercept_lin <- coef(linear_mod)[1]
fast_food_spend_lin <- coef(linear_mod)[2]
coffee_drinker_lin <- coef(linear_mod)[3]

# --- 2) Extract coefficients from logit regression model
intercept_logit <- coef(logit_mod)[1]
fast_food_spend_logit <- coef(logit_mod)[2]
coffee_drinker_logit <- coef(logit_mod)[3]

# --- Estimate predicted data from different models   
heart_data <- heart_data %>%
  mutate(
    # Convert outcome variable to factor
    heart_disease_factor = factor(heart_disease, 
                                  labels = c("No Disease (Y=0)", "Disease (Y=1)")),
    # 1) Linear model prediction
    lin_pred = intercept_lin + fast_food_spend_lin * fast_food_spend, 
    # 2) Logit model prediction
    logit_pred = intercept_logit + fast_food_spend_logit * fast_food_spend,  coffee_drinker,
    # 3) Convert logit to probability (logistic model prediction)
    logistic_pred = 1 / (1 + exp(-logit_pred)) 
    ) %>%
  arrange(fast_food_spend)   

• Besides the poor fit, recall linear regression models implies certain assumptions:
  • Linear relationship between Y and the predictors X1 and X2
    
  • Residuals must be 1) approximately normally distributed and 2) uncorrelated
  • Homoscedasticity: residuals should have constant variance
  • Non collinearity: predictor variables should not be highly correlated with each other

• Assumptions that are not met in this case

# diagnostic plots
par(mfrow=c(2,2))
plot(linear_mod)

• Table 1 (next) shows the model output, with the coefficient estimate for each predictor.
  • The broom::tidy function converts the model summary into a more readable data frame.
  • The odds ratio (= exponentiated coefficient estimate) is more interpretable than the coefficient itself.

# Convert model's output summary into data frame
heart_model_coef <- broom::tidy(heart_model) %>% 
  # improve readability of significance levels
  dplyr::mutate('signif. lev.' = case_when(
    `p.value` < 0.001 ~ "***",
    `p.value` < 0.01 ~ "**",
    `p.value` < 0.05 ~ "*",
    TRUE ~ ""))%>%
  # add odds ratio column
  dplyr::mutate(odds_ratio = exp(estimate)) %>%
  dplyr::relocate(odds_ratio, .after = estimate) %>%
  dplyr::mutate(across(where(is.numeric), ~ round(.x, 4))) %>%
  # format as table
  knitr::kable() %>% 
  # reduce font size
  kable_styling(font_size = 20) %>% 
  # add table title
  kableExtra::add_header_above(c("Logistic Regression Analysis of Heart Disease Risk Factors"= 7))

heart_model_coef

• Intercept: gives the log-odds of heart disease when all predictor variables are zero. This is not generally interpreted, but the highly negative value suggests very low probability of heart disease in the sample of reference.
  • (If interpreted) the intercept term -11.0554 would translate as probability \(P = e^{-11.0554} / (1 + e^{-11.0554}) = 0.00002\) or \(0.002\%\).
  • …which means: when \(\text{coffee_drinker} = NO\), and \(\text{fast_food_spend} = 0\), and \(\text{income} = 0\), the probability of heart disease is as low as \(0.002\%\).

• Income: Based on a p-value = 0.8182, we conclude that income is not significantly associated with heart disease. (Anyhow, the odds ratio of ≈1 would suggest no change in odds based on income.)

• The positive coefficient of Fast Food Annual Spending (US$), 0.0024 (highly statistically significant as p-value = 0.0000), suggests a positive association with heart disease.

• This means that for each additional \(\Delta X_1 = +1 \; US\$\) spent on fast food annually:

  • the log(odds) of heart disease increases by: \(\Delta \log(\text{odds}) = \beta_{1} \times \Delta X_1 = 0.0024 \times 1 = 0.0024\)

  • the odds of heart disease (with spending) increase by a factor of: \(OR = e^{0.0024} \approx 1.0024\) (compared to without spending).

• The probability of heart disease is computed using the logistic function: \(P_{HD=1} = \frac{e^{\beta_0 + (\beta_1 \times X_1)}}{1 + e^{\beta_0 + (\beta_1 \times X_1)}}\)
  • Solving for \(\beta_0 = -11.0554\), \(\beta_1 = 0.0024\) and \(X_1 = 1\) gives:
  \[P_{HD=1} = \frac{e^{-11.0554 + (0.0024 \times 1)}}{1 + e^{-11.0554 + (0.0024 \times 1)}} = 0.00001583981\]

which indicates a probability of heart disease is as low as \(0.00159\%\) when fast_food_spend \(= 1\; US\$\).

• However, considering that this is annual spending, we should use a more adequate scale!

• For example, for an additional \(\Delta X_1 = +100 \; US\$\) spent on fast food annually:

  • the log(odds) of heart disease increases by: \(\Delta \log(\text{odds}) = \beta_{1} \times \Delta X_1 = 0.0024 \times 100 = 0.24\)

  • the odds of heart disease (with spending) increase by a factor of: \(OR = e^{0.24} \approx 1.271\) (compared to without spending).

• The probability of heart disease is computed using the logistic function: \(P_{HD=1} = \frac{e^{\beta_0 + (\beta_1 \times X_1)}}{1 + e^{\beta_0 + (\beta_1 \times X_1)}}\)

• Solving for \(\beta_0 = -11.0554\), \(\beta_1 = 0.0024\) and and \(X_1 = 100\) gives:
  \[P_{HD=1} = \frac{e^{-11.0554 + (0.0024 \times 100)}}{1 + e^{-11.0554 + (0.0024 \times 100)}} \approx 0.00002008\]

  which (this time) indicates a (still very low probability) of heart disease at \(0.002008\%\) when fast_food_spend \(= 100\; US\$\).
  

• The negative coefficient of Coffee Drinker (=YES), -0.7296 expressed in log-odds, means that coffee drinking is associated with lower odds of having heart disease.

• Transforming the estimate into odds ratio, we obtain \(0.48 = 48\%\), which tells that coffee-drinkers have 0.48 (or 48%) lower odds than non-coffee drinkers to experience heart disease (holding all other predictors fixed). \[ \text{Odds Ratio} = e^{-0.7296} = 0.48 \]

  • Alternatively, we could say that non-coffee drinkers have \(1-0.48 = 0.52 = 52\%\) higher odds of having heart disease compared to coffee drinkers.

• This effect is statistically significant as p-value = 0.0122.

Name: Biopsy Data on Breast Cancer Patients
Documentation: See reference on the data downloaded and conditioned for R here https://cran.r-project.org/web/packages/MASS/MASS.pdf
Sampling details: This breast cancer database was obtained from the University of Wisconsin Hospitals, Madison from Dr. William H. Wolberg. He assessed biopsies of breast tumours for 699 patients up to 15 July 1992; each of nine attributes has been scored on a scale of 1 to 10, and the outcome is also known. The dataset contains the original Wisconsin breast cancer data with 699 observations on 11 variables.

The biopsy data contains 699 observations of 9 continuous variables, V1, V2, …, V9.

The dataset also contains a character variable: id, and a factor variable: class, with two levels (“benign” and “malignant”).

# check variable types
str(biopsy)

'data.frame':   699 obs. of  11 variables:
 $ ID   : chr  "1000025" "1002945" "1015425" "1016277" ...
 $ V1   : int  5 5 3 6 4 8 1 2 2 4 ...
 $ V2   : int  1 4 1 8 1 10 1 1 1 2 ...
 $ V3   : int  1 4 1 8 1 10 1 2 1 1 ...
 $ V4   : int  1 5 1 1 3 8 1 1 1 1 ...
 $ V5   : int  2 7 2 3 2 7 2 2 2 2 ...
 $ V6   : int  1 10 2 4 1 10 10 1 1 1 ...
 $ V7   : int  3 3 3 3 3 9 3 3 1 2 ...
 $ V8   : int  1 2 1 7 1 7 1 1 1 1 ...
 $ V9   : int  1 1 1 1 1 1 1 1 5 1 ...
 $ class: Factor w/ 2 levels "benign","malignant": 1 1 1 1 1 2 1 1 1 1 ...

There is one incomplete variable V6 = “bare nuclei” with 16 missing values.

• remember the package skimr for exploring a dataframe?

# check if vars have missing values
biopsy %>% 
  # select only variables starting with "V"
  skimr::skim(starts_with("V")) %>%
  dplyr::select(skim_variable, 
                n_missing)

# A tibble: 9 × 2
  skim_variable n_missing
  <chr>             <int>
1 V1                    0
2 V2                    0
3 V3                    0
4 V4                    0
5 V5                    0
6 V6                   16
7 V7                    0
8 V8                    0
9 V9                    0

We can decide what to do in these cases (informed by our knowledge of the dataset):

• Option 1) We drop the observation with incomplete data (i.e. with missing values for V6 = “bare nuclei”) with 16 missing values.

# remove rows with missing values
biopsy_drop <- biopsy %>% 
  dplyr::filter(!is.na(V6))

mean(biopsy_drop$V6)

[1] 3.544656

• Option 2) We impute the missing values with the mean of the variable V6 = “bare nuclei”.

# impute missing values with the median of the variable
biopsy_impute <- biopsy %>% 
  dplyr::mutate(V6 = ifelse(is.na(V6), median(V6, na.rm = TRUE), V6))

mean(biopsy_impute$V6)

[1] 3.486409

• Biopsied cells of 700 breast cancer tumors, used to determine if the tumors were benign or malignant.
• This determination was based on 9 characteristics of the cells, ranked from 1(benign) to 10(malignant):
  • 1) Clump Thickness – How the cells aggregate. If monolayered they are benign and if clumped on top of each other they are malignant
  • 2) Uniform Size – All cells of the same type should be the same size.
  • 3) Uniform Shape If cells vary in cell shape they could be malignant
  • 4) Marginal Adhesion – Healthy cells have a strong ability to stick together whereas cancerous cells do not
  • 5) Single Epithelial Size – If epithelial cells are not equal in size, it could be a sign of cancer
  • 6) Bare nuclei – If the nucleus of the cell is not surrounded by cytoplasm, the cell could be malignant
  • 7) Bland Chromatin – If the chromatin’s texture is coarse the cell could be malignant
  • 8) Normal Nucleoli – In a healthy cell the nucleoli is small and hard detect via imagery. Enlarged nucleoli could be a sign of cancer
  • 9) Mitosis – cells that multiply at an uncontrollable rate could be malignant

• The explanatory variable(s) (clump_thickness, …, mitosis) can be renamed for better readability
• The observations with missing values (bare_nuclei) are removed for simplicity
  
• Patient ID (id) can be dropped as it is not used in this analysis

# Create a clean version of the dataset
biopsy_clean <- biopsy %>%
  # rename the columns (new = old)
  rename(
    id                = ID,
    clump_thickness   =  V1,
    uniform_size      =  V2,
    uniform_shape     =  V3,
    marginal_adhesion =  V4,
    single_epith_size =  V5,
    bare_nuclei       =  V6,
    bland_chromatin   =  V7,
    normal_nuclei     =  V8,
    mitosis           =  V9,
    class             =  class) %>% 
  # remove rows with missing values
  na.omit(bare_nuclei) %>% 
  # remove the id column
  select(-id)

# check the structure of the dataset
paint::paint(biopsy_clean)

data.frame [683, 10]
clump_thickness   int 5 5 3 6 4 8
uniform_size      int 1 4 1 8 1 10
uniform_shape     int 1 4 1 8 1 10
marginal_adhesion int 1 5 1 1 3 8
single_epith_size int 2 7 2 3 2 7
bare_nuclei       int 1 10 2 4 1 10
bland_chromatin   int 3 3 3 3 3 9
normal_nuclei     int 1 2 1 7 1 7
mitosis           int 1 1 1 1 1 1
class             fct benign benign benign benign benign ma~

• We fit a logistic regression model to the biopsy_train data using:
  • the glm function with argument family = binomial to specify the logistic regression model;
  • and with Class ~ . to specify an initial model that uses all the variables as predictors (backward elimination approach).

# Building initial model 
model = stats::glm(class ~ . , family = binomial, data=biopsy_train)

• Table 2 shows the model summary, with the coefficient estimate for each predictor.
  • the broom::tidy function converts the model summary into a data frame.

broom::tidy(model) %>% 
  mutate('Sign.lev' = case_when(
    `p.value` < 0.001 ~ "***",
    `p.value` < 0.01 ~ "**",
    `p.value` < 0.05 ~ "*",
    TRUE ~ ""))%>%
  mutate(across(where(is.numeric), ~ round(.x, 4))) %>%
  knitr::kable() 

In logistic regression, the coefficients are in the form of the natural logarithm of the odds of the response event happening (i.e. \(Y_i = 1\)):

\[logit(p_i) = \ln\left( \frac{p_i}{1-p_i} \right) = -9.5063 + 0.3935 \times Clump\_Thickness + ... + 0.5065 \times Mitosis\]

However, with some algebraic transformation, the logit function can be inverted to obtain the probability of the response event happening as a function of the predictors:

\[p_i = \frac{1}{1 + e^{-(-9.5063 + 0.3935 \times Clump\_Thickness + ... + 0.5065 \times Mitosis)}}\]

This equation represents the logistic regression model’s best-fit line.

• The model includes all variables, but we could make it more parsimonious by removing variables that are not significant!
• We use a statistic called the Akaike Information Criterion (AIC) to compare models.
  • AIC’s calculation gives a penalty for including additional variables.
  • The model with the lowest AIC is considered the best.

# For example let's fit a model without the variable `uniform_size`
model2 = glm(class~ .-uniform_size, family = binomial, data=biopsy_train)

According to the AIC values, the model2 seems better (AIC is lower).

# Compare the AIC values of the 2 models
tibble(Model = c("model", "model2"), 
       AIC = c(AIC(model), AIC(model2) )) %>% 
  kable()

• The MASS package’s function stepAIC enables to perform a systematic model selection (by AIC):
  • The direction argument specifies the direction of the stepwise search.
  • The trace argument (if set to TRUE) prints out all the steps.
• The best_model has removed these variables:
  • uniform_shape
  • single_epith_size
  • normal_nuclei
• The best_model has the lowest AIC value (from 100 to 98.5), despite a higher Residual Deviance than the full model (from 80 to 80.6), albeit by a very slight amount.

# Select the best model based on AIC
best_model <- MASS::stepAIC(model, direction = "both", trace = FALSE)

# Compare the AIC values of full and best model
tibble(Model = c("model", "best_model"), 
       AIC = c(AIC(model), AIC(model2)),
       Deviance = c(deviance(model), deviance(best_model))) %>% kable()

A 3D scatterplot of observations shows the first 3 principal components’ scores.

• For this one, we need the scatterplot3d() function of the scatterplot3d package;
• The color argument assigned to the Label variable;
• To add a legend, we use the legend() function and specify its coordinates via the xyz.convert() function.

# 3D scatterplot ...
plot_3d <- with(PC_scores, 
                scatterplot3d::scatterplot3d(PC_scores$PC1, 
                                             PC_scores$PC2, 
                                             PC_scores$PC3, 
                                             color = as.numeric(Label), 
                                             pch = 19, 
                                             main ="Figure 3: 3D Scatter Plot", 
                                             xlab="PC1",
                                             ylab="PC2",
                                             zlab="PC3"))

# ... + legend
legend(plot_3d$xyz.convert(0.5, 0.7, 0.5), 
       pch = 19, 
       yjust=-0.6,
       xjust=-0.9,
       legend = levels(PC_scores$Label), 
       col = seq_along(levels(PC_scores$Label)))

Biplots have two key elements: scores (the 2 axes) and loadings (the vectors). As in the scores plot, each point represents an observation projected in the space of principal components where:

• Biopsies of the same class are located closer to each other, which indicates that they have similar scores referred to the 2 main principal components;
• The loading vectors show strength and direction of association of original variables with new PC variables.

As expected from PCA, the single PC1 accounts for variance in almost all original variables, while V9 has the major projection along PC2.